
program define qte, eclass sortpreserve
	version 7.0

	syntax varlist(min=2 max=2) [if] [in] [aweight fweight /] , Quantiles(numlist) [ model(string) covariates(varlist numeric) nodid nodit nodic firpo(string asis) PSCOREname(varlist min=1 max=1) treatmnt(varlist min=1 max=1) INFlatecov(string asis) noKd dots noJcumul Generate(string) noTG noCG] 

	if ("`tg'`cg'"=="notgnocg") {
		di _new " -> You can't use both the -notg- and -nocg- options at once." 
		exit
	} 
	else if "`tg'"=="notg" {
		local dit "nodit"
		local did "nodid"
	}
	else if "`cg'"=="nocg" {
		local dic "nodic"
		local did "nodid"
	}

	if ("`dit'`dic'`did'"=="noditnodicnodid") {
		di _new " -> You can't use all three of nodit, nodic, and nodid at once." 
		exit
	} 

	local cumul
	if "`jcumul'"=="nojcumul" {
		local cumul "cumul"
	}
	else if "`jcumul'"=="" {
		local cumul "jcumul"
	}

	*forcing to use cumul, since jcumul code not working right now
	*local cumul "cumul"

	if "`inflatecov'"=="" { 
		local inflatecov 1
	}

	if "`model'"=="" & "`covariates'"~="" { /* default model */
		local model "logit"
	} /* end of default model */

	/* giving weight var a mnemonic name */
	if "`weight'" ~= "" { 
		local weightexp "`exp'" 
	}

	/* mark sample */
	marksample touse 
	qui count if `touse'
	local Nboth = _result(1)

	/* parse varlist */
	local depvar : word 1 of `varlist'
	local treat  : word 2 of `varlist'

	/* check that treatment dummy is in fact binary */
	qui count if `treat' ~= 0 & `treat' ~= 1 & `touse'
	if _result(1) ~= 0 { /* exit bc treatment var is not binary in sample */
		di in red "Treatment variable '`treat'' is not binary."
		error 
	} /* end of exit bc treatment var is not binary in sample */
 
	/* run firpo if asked for it */
	tempvar pscore pscorewt
	if "`firpo'"~="" {
		`firpo' /* runs the model, which has to be fully specified by the user. compound quotes bc of asis option */
		qui predict double `pscore' if `touse'
		qui gen double `pscorewt' = `treatmnt'/`pscore' + (1-`treatmnt')/(1-`pscore') if `touse'
		if "`pscorename'"~="" { qui gen double `pscorename' = `pscore' if `touse' }
		su `pscorewt'
	} 
	else { 
		*user didn't ask for firpo model, so we define pscorewt var as 1 for every obs
		qui gen `pscorewt' = 1 if `touse'
	}

	/* this is the weight we actually use */
	/* MB changes this on 8/11 to use firpo with final weight if firpo called */
	tempvar finalweight
	if "`weight'" ~= "" { 
		qui gen double `finalweight' = `weightexp' * `pscorewt' if `touse' 
	}
	else if "`firpo'" != ""	    { /* firpo weights */ 
		qui gen double `finalweight' = `pscorewt'      if `touse' 
	}
	else 			    { /* constant weights */ 
		qui gen double `finalweight' = 1               if `touse' 
	}


	* 8/11 debug mb
	*di "finalweight"
	*su `finalweight'


	tempvar cumvar
	qui gen double `cumvar' = . /* initialization */

	*creating ecdf var => `cumvar'
	/* used to do cumul using bysort `treat'
	   but this is no good: need to do it explicitly or stata screws up ties within bysortvar when we use cumul with bysort
	   prev code was ok even so, b/c of the sort (below) on `treat' `cumvar' in step before getting quantiles
	   fixed it anyway to stay clean, and in case that code ever changes
	*/
	di _new "Estimating CDF (using `cumul' command)..."

	*use jcumul, and hence N=n0+n1 in denom
	if "`cumul'"=="jcumul" {

*NEED TO ADD CODE TO IGNORE GROUP WHEN CALL WITH notg OR nocg option

		tempvar dummy
                qui jcumul `depvar' if `touse' `in' [aw=`finalweight'], gen(`dummy') by(`treat')
		qui replace `cumvar' = `dummy'
		qui drop `dummy'
	}
	else {
		forvalues t = 0/1 { /* looping over treatment dummy to do cumul thing */

			if ("`t'`cg'"=="0nocg") {
				continue
			}
			else if ("`t'`tg'"=="1notg") {
				continue
			}

			tempvar dummy`t'
			qui cumul `depvar' if `touse' & `treat'==`t' `in' [aw=`finalweight'], gen(`dummy`t'')
			qui replace `cumvar' = `dummy`t'' if `treat'==`t'
			qui drop `dummy`t''
		} /* end of looping over treatment dummy to do cumul thing */
	}
	*end of else condn on which cumul command to use

	local M : word count `quantiles'	/* number of quantiles estimated */
	local M = `M' + 1 /* add 1 for the means */

	tempname b0 b1 f0 f1 v0 v1
	mat `b0' = J(1,`M',0)
	mat `b1' = J(1,`M',0)
	mat `f0' = J(1,`M',0)
	mat `f1' = J(1,`M',0)
	mat `v0' = J(`M',`M',0)
	mat `v1' = J(`M',`M',0)


	*setting up vars used in kdensity thing
	if "`kd'"=="" {
		tempvar kd0 kd1 kdatvar0 kdatvar1 
		qui for new `kdatvar0' `kdatvar1' : qui gen double X  = .
	} /* end of setting up vars used in kdensity thing */

	forvalues t = 0/1 { /* looping over treatment dummy */

		if ("`t'`cg'"=="0nocg") {
			continue
		}
		else if ("`t'`tg'"=="1notg") {
			continue
		}

		/* counter for quantile we're on */
		local n 0
		foreach q in `quantiles' { /* looping over quantiles to run qregs */

			local n = `n' + 1 /* incrementing quantile counter */

			/* making local macro to store value of each quantile (`q`M'') */
			local q`n' = `q'

			if `q' <= 0 | `q' >= 100 {
				di in red "quantiles(`q') out of range"
				exit 198
			}

			/* using cumul approach--should bstrap the cov matrix, tho we will provide kdensity-based cov estimates */
			qui sort `treat' `cumvar'
			qui sum `depvar' `in' [aw=`finalweight'] if /*
				*/ (`touse'==1) & (`cumvar' >=(`q'/100)) & (`cumvar'[_n-1]<(`q'/100)) & (`treat'==`t')
		        qui mat `b`t''[1,`n'] = _result(3)	/* estimated quantile 		*/
			qui count `in' if `touse' & `treat'==`t'
			local numobs`t' = _result(1)

			/* making names */
			if `t'==0 { /* control group names */
				local namesc "`namesc' c:q`q'"
			}
			else {
				local namest "`namest' t:q`q'"
				local namesd "`namesd' d:q`q'"
			} /* end of names */

			*preparing for density business below
			if "`kd'"=="" {
				qui replace `kdatvar`t'' = `b`t''[1,`n'] if _n== `n'
			} /* end of preparing for density business below */

			if ("`dots'"~="") & (`n'/2==int(`n'/2)) { /* do every other b/c treatment loop is outside quantiles loop */
				di "." _continue
			}
		} /* end of looping over quantiles */

		*getting density estimates
		if "`kd'"=="" {
			kdensity `depvar' [aw=`finalweight'] if `touse' & `treat'==`t' , gen(`kd`t'') at(`kdatvar`t'')

			local mm1 = `M'-1 /* last one is the mean, remember */
			forvalues r = 1/`mm1' { /* looping over row quantiles to get density estimates and fix cov matrix */
				forvalues c = 1/`mm1' { /* looping over col quantiles to fix cov matrix */
	
					*getting density estimates
					local rdensity = `kd`t''[`r']		/* getting approp kdensity estimate for row term */
					local cdensity = `kd`t''[`c']		/* getting approp kdensity estimate for col term */
	
					*getting part of cov term assoc w/the particular quantiles
					*(see thm 3.1 of k & b (1982) ema paper)
					local qpart = min(`q`r'', `q`c'')/100 - `q`r''*`q`c''/10000
	
			                mat `v`t''[`r',`c'] = ( `qpart' / (`rdensity'*`cdensity') ) * `inflatecov'
			                mat `v`t''[`c',`r'] = ( `qpart' / (`rdensity'*`cdensity') ) * `inflatecov'
	
				} /* end of looping over col quantiles to get density estimates and fix cov matrix */
			} /* end of looping over row quantiles to get density estimates and fix cov matrix */	

			mat `v`t'' = `v`t''*(1/`numobs`t'')
		} /* end of getting density estimates */

		/* now do means */
		qui ci `depvar' if `touse' & `treat'==`t' [aw=`finalweight']
		mat `b`t''[1,`M']	= $S_3
		mat `v`t''[`M',`M'] 	= ($S_4)^2

		if `t'==0 { /* control group names for mean */
			local namesc "`namesc' c:mean"
		}
		else {
			local namest "`namest' t:mean"
			local namesd "`namesd' d:mean"
		} /* end of names */

		/* naming matrices -- have to do this outside quantile loop */
		if `t'==0 { /* control group names */
			mat colnames `b`t'' = `namesc'
			mat rownames `v`t'' = `namesc'
			mat colnames `v`t'' = `namesc'
		}
		else {
			mat colnames `b`t'' = `namest'
			mat rownames `v`t'' = `namest'
			mat colnames `v`t'' = `namest'
		} /* end of names */
	} /* end of looping over treatment dummy */ 


	/* now doing the wrap-up */
	tempname b v bd vd bdpost vdpost zero

	/* dif, i.e. qte */
	mat `bd'  = `b1' - `b0'
	mat colnames `bd' = `namesd'

 	mat `vd' = `v1' + `v0'
 	mat rownames `vd' = `namesd'
	mat colnames `vd' = `namesd'

	*use this mat for adding in blocks in constructing full cov matrix (with control, treatment, and dif)
	mat `zero' = J(`M',`M',0)	

	di
	di in gr "Quantile treatment effects" /*
			*/ _col(54) "Number of T obs =" in ye %8.0g `numobs1'

	if "`dit'`dic'" == "" { /* want tg and cg results */

		if "`did'"=="" {
			mat `b' = `b1' , `b0', `bd'
			mat `v' = (`v1' , `zero' , `v1') \ (`zero' , `v0', -1*`v0') \ (`v1', -1*`v0', `vd')
			mat rownames `v' = `namest' `namesc' `namesd'
			mat colnames `v' = `namest' `namesc' `namesd'
			di "(Treated and control quantiles and QTE)" _col(54) /*
				*/	    "Number of C obs =" in ye %8.0g `numobs0'
		} /* end `did'=="" */
		else {
			mat `b' = `b1' , `b0'
			mat `v' = (`v1' , `zero') \ (`zero' , `v0')
			mat rownames `v' = `namest' `namesc' 
			mat colnames `v' = `namest' `namesc' 
			di "(Treated and control quantiles)" _col(54) /*
				*/	    "Number of C obs =" in ye %8.0g `numobs0'
		} /* end `did'=="" */
	}
	else if "`dic'" == "" { /* want only the control and dif results */

		if "`did'"=="" {
			mat `b' = `b0', `bd'
			mat `v' = (`v0', -1*`v0') \ ( -`v0', `vd')
			mat rownames `v' = `namesc' `namesd'
			mat colnames `v' = `namesc' `namesd'
			di "(Control quantiles and QTE only)" _col(54) /*
				*/	    "Number of C obs =" in ye %8.0g `numobs0'
		} /* end `did'=="" */
		else {
			mat `b' = `b0'
			mat `v' = (`v0') 
			mat rownames `v' = `namesc' 
			mat colnames `v' = `namesc' 
			di "(Control quantiles and QTE only)" _col(54) /*
				*/	    "Number of C obs =" in ye %8.0g `numobs0'
		} /* end `did'=="" */
	}
	else if "`dit'" == "" { /* want only the treatment and dif results */

		if "`did'"=="" { 
			mat `b' = `b1', `bd'
			mat `v' = (`v1', `v1') \ ( `v1', `vd')
			mat rownames `v' = `namest' `namesd'
			mat colnames `v' = `namest' `namesd'
			di "(Treated quantiles and QTE only)" _col(54) /*
				*/	    "Number of C obs =" in ye %8.0g `numobs0'
		} /* end `did'=="" */
		else {
			mat `b' = `b1'
			mat `v' = (`v1')
			mat rownames `v' = `namest' 
			mat colnames `v' = `namest' 
			di "(Treated quantiles and QTE only)" _col(54) /*
				*/	    "Number of C obs =" in ye %8.0g `numobs0'
		} /* end `did'=="" */
	}
	else if ("`dit'"~="" & "`dic'"~="") { /* want only the dif results */
		mat `b' = `bd'
		mat `v' = `vd'

		/* losing the equation name */
		local names : rownames `v'
		mat colnames `b' = `names'
		mat colnames `b' = _:
		mat rownames `v' = `names'
		mat rownames `v' = _:
		mat colnames `v' = `names'
		mat colnames `v' = _:

		di "(QTE only)" _col(54) /*
			*/	    "Number of C obs =" in ye %8.0g `numobs0'
	} /* end of subsets of the results */


        estimates post `b' `v' , esample(`touse') depname(`depvar') 

	/* quantile estimates */
	estimates matrix b0 `b0'
	estimates matrix b1 `b1'
	estimates matrix bd `bd'

	/* variance estimates */
	estimates matrix v0 `v0'
	estimates matrix v1 `v1'
	estimates matrix vd `vd'

	/* scalars */
	if "`cg'"=="" {
		estimates scalar N0 = `numobs0'
	}

	if "`tg'"=="" {
		estimates scalar N1 = `numobs1'
	}

	/* storing cdf estimates if requested */
	if "`generate'"~="" { 
		qui rename `cumvar' `generate'
	}

	di
	estimates display

end

